Phase coexistence of cluster crystals: beyond the Gibbs phase rule 
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We report a study of the phase behavior of multiple-occupancy crystals through simulation. We 
argue that in order to reproduce the equilibrium behavior of such crystals it is essential to treat 
the number of lattice sites as a constraining thermodynamic variable. The resulting free-energy 
calculations thus differ considerably from schemes used for single-occupancy lattices. Using our 
approach, we obtain the phase diagram and the bulk modulus for a generalized exponential model 
that forms cluster crystals at high densities. We compare the simulation results with existing 
theoretical predictions. We also identify two types of density fluctuations that can lead to two 
sound modes and evaluate the corresponding elastic constants. 

PACS numbers: 64.70.Kb, 64.70.Dd, 82.30.Nr, 83.80.Rs 



At finite temperatures all crystals contain point de- 
fects. This means that the ratio between N, the number 
of particles, and N c , the number of unit cells is not fixed 
by geometry. In the language of Ref. [l[ , N c is a "con- 
strained" thermodynamic variable. A general variation 
of the Hclmholtz free energy for a one-component crystal 
can be written as: 



dF = -SdT - PdV + ndN + fi c dN c 



(1) 



where S is the entropy, T the absolute temperature, P 
the pressure, V the volume, /x the chemical potential of 
the constituent particles, and \x c the "cell chemical po- 
tential" conjugate to the number of unit cells. If N c is 
free to change, it will take on a value such that fi c = 0, 
hence its value is a function of N, V, and T. In simple 
crystals, the equilibrium concentration of point defects is 
usually so low that their effect on the phase behavior is 
negligible @, For instance, at melting the chemical 
potential of a hard-sphere crystal with vacancies roughly 
differs by as little as 10 -3 £:bT from that of a defect-free 
crystal for which N c = N Q. 

Interestingly, the situation is dramatically different in 
"cluster crystals" @, i, i, 0, i, 0. These unusual crys- 
talline materials can have a number of particles per lat- 
tice site much larger than one. Such solids form in 
systems of particles that interact via a bounded, short- 
ranged and purely repulsive pair potential whose Fourier 
transform has negative regions. The effect of fi c on the 
phase behavior then becomes all important, which has 
profound consequences for the numerical study of their 
phase transitions. The reason is that in almost all simu- 
lations involving crystals, the average number of particles 
per unit cell is fixed at the outset of the simulation. Af- 
ter that, a change in the density p = N/V of the system 
may still change P and fi but, as the ratio N/N c is fixed, 
He will in general not be zero. Hence, conventional sim- 
ulations do not probe the lowest free-energy state of the 



crystal. At constant P and T, a small variation in Gibbs 
free energy G = F + PV is of the form fidN + p, c dN c . 
If we fix the ratio n c = N/N c , then both /1 and fi c are 
constant, so we can integrate to obtain: 



and hence 



G = Nfi + N c n c 



N c fi c = F + PV - /J.N. 



(2) 



(3) 



For a given N, V, T, and N C) we can use Monte Carlo 
(MC) simulations to compute F, P, and /z As all 

quantities on the right-hand side of Eq. [3] can be de- 
termined numerically, whilst N c is known, we can also 
compute fi c . This is important because the condition 
for phase coexistence involving cluster crystals requires 
equality of \i, P, and T in the coexisting phases and of 
|U C = in all crystalline phases. This latter condition is 
not normally considered in the discussion of the Gibbs 
phase rule. However, in his original formulation, Gibbs 
does allow for the possible existence of other thermody- 
namic "fields" " 



n|. 



in addition to /1, P and T 
As an application of this approach we consider the nu- 
merical determination of the phase diagram for the gener- 
alized exponential model (GEM-n) = ee^^ rij ^ a ^" 
with n = 4, where e and a determine respectively the en- 
ergy and the length scales. For convenience, we set them 
to unity and consider only reduced units from this point 
forward. For n > 2, this system is known to form cluster 
solids at high densities 0, H, H, HH . Its phase diagram is 
known qualitatively, but not quantitatively: at high T, 
the fluid first freezes into a multiply-occupied BCC phase 
that transforms into a multiply-occupied FCC phase at 
higher densities. By contrast, upon compression at low 
T, the system undergoes a "normal" freezing transition 
to a single-occupancy FCC crystal; clustering only sets 
in upon further compression of the solid. 
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FIG. 1: Comparison of the DFT [5|] (dashed lines) and 
SAMC 7] (dash-dotted lines) T-p phase diagrams with the 
simulation results (points) for the GEM-4 model. The gray 
zone highlights the phase coexistence region. Upper inset: P- 
T simulation phase diagram (points) with the corresponding 
Clausius-Clapeyron tangents to the coexistence curve (black 
lines). The solid curves are guides for the eye. Lower in- 
set: shifted free energy curves F = F/VT — 26p (to enhance 
visibility) at T = 0.2 for the liquid (dash-dotted gray line), 
BCC crystal (solid gray line), and FCC crystal (dashed black 
line), along with the common tangent construction (solid 
black lines) and the coexistence densities (dots and drop-down 
lines). 



We perform discretized-space constant- NVT MC simu- 
lations 0, EH for 2000-5000 particles in the temperature 
regime where multiple occupancy of the crystal lattice 
sites is expected. We determine the value of p c = N c /V 
such that (x c = for every (p, T) point by starting the 
fixed N c simulations with a reasonable guess for N and V 
and iterating until the correct values of N and V are lo- 
cated. Via the common tangent construction, we obtain 
the coexistence densities from the resulting free energy 
curves. We note that, because the number of particles 
per lattice site is free to fluctuate, we cannot use the 
Einstein-crystal method to compute the free energy of 
the solid 10, 14]. Rather, we perform a thermodynamic 
integration from a reference state of ideal-gas particles 
that move in potential wells centered around the lattice 
We also note that for particles that form cluster 



15] 



sites 

solids, the Widom particle- insertion method provides an 
efficient tool to determine the chemical potential, even in 
the dense solid 0, [l6| . One could even think of perform- 
ing a kind of Gibbs-ensemble simulation where two sys- 
tems exchange both particles and volume 17j . However, 
such a simulation would not locate the correct coexistence 
point, precisely because the Gibbs-ensemble method does 
not impose the condition p c = in all solid phases. That 
is why we have to follow the rather elaborate route via 
Eq. [3] to locate the points where the different phases co- 
exist. 



In Fig. Q] we compare the T-p phase diagram for the 
GEM-4 model with the corresponding density-functional 
theory (DFT) predictions of Ref. @]. In the same figure, 
we also show the estimate of the freezing transition based 
on the results of simulated annealing MC (SAMC) 0,(1]. 
For the liquid-BCC transition, the liquidus line predicted 
by DFT is indistinguishable from the "exact" simulation 
results and the solid line is only slightly off. The SAMC 
results, although qualitatively correct, predict a liquid- 
BCC density gap that is too narrow. This is probably 
due to finite-size artefacts in Refs. 0, Q ■ A crucial test 
of the accuracy of DFT is the prediction of the location 
of the BCC-FCC transition. As the densities and free 
energies of these two phases are very close (Fig. [T] lower 
inset), minor inaccuracies in the theory should have a no- 
ticeable effect on the prediction of the transition point. 
Indeed, we find that even though the DFT free-energy 
predictions are only off by a small amount (not shown) , 
the location of the phase transition is shifted by roughly 
10% in p. The P-T phase diagram is shown in the up- 
per inset, where the various state points are accompanied 
by tangents to the coexistence curves obtained from the 
simulation free and internal energy results by use of the 
Clausius-Clapeyron relation [103]. Inspection of the di- 
agram suggests a liquid-BCC-FCC triple point around 
T t ~ 0.15, which is much lower than the DFT prediction 
of T t s» 0.4 The dramatic shift follows from the small 
difference between the slopes of the liquid-BCC and the 
BCC-FCC coexistence curves, leaving the location of the 
triple point very sensitive to any modification of the lat- 
ter. 

It is interesting to understand the reason for the failure 
of DFT to predict the location of the solid-solid transi- 
tion. One of the core consequences of the DFT treatment 
of Ref. H is that the volume of the unit cell in a cluster 
crystal v c = V/N c is independent of density, so that the 
average site occupation n c cx p. This feature is known to 
break down in low-density crystals and is also found 
here to be slightly inaccurate at intermediate tempera- 
tures and densities. Though the linear relationship holds, 
the proportionality is shifted by a constant. As can be 
gathered from the inset of Fig. [2j this leads to a non-zero 
value of (dv c /dp)T for equilibrium states, though the ef- 
fect vanishes with increasing T and p. This suggests that 
the DFT approximation is asymptotically valid. At in- 
termediate densities, this correction though small might 
be sufficient to explain the discrepancy between the DFT 
and the numerical results. 

Since thermodynamic equilibrium is only obtained 
when p c = 0, it would appear that once the equilibrium 
points are found, all references to the artificial p c field can 
be disregarded. Yet, for quantities involving the second 
derivative of the constrained free energy, such as the bulk 
modulus B = V (§p§J , it cannot be neglected un- 
less one already has the complete equilibrium free energy 
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FIG. 2: Bulk modulus results from direct differentiation of 
the free energy for three different temperatures in the stable 
crystal structures [BCC (gray) and FCC (black) at T = 0.5 
(solid line), T = 0.8 (dashed lined), and T = 1.1 (dash-dotted 
line)], along with the values at three state points (stars) com- 
puted using Eq. [3] The virial contribution to B is shown for 
reference (crosses). The breakdown of B is also given in Ta- 
ble U Inset: variation of the lattice volume with density in 
equilibrium for these same systems. DFT presents this quan- 
tity as zero. 



curve at hand. In simulations of single-occupancy crys- 
tals, -B v i r = —V(§£-) NTN can be computed directly for 
a given state point throu gh an approach similar to the 
virial calculation of P [la , [l9| . For cluster crystals, how- 
ever, the artificial system conditions further modify the 
bulk modulus as 



B = B v 



EL 

n r 



P_ 



dv c 
dp 



, (4) 



where the partial derivatives are evaluated around an 
equilibrium state point. The virial contribution corre- 
sponds to a quenched system where particle rearrange- 
ments are not possible, so it is an upper bound to 
B = B vir — B COII . The results for different state points 
are compared in Fig. [3] to the values obtained by direct 
numerical differentiation of the equilibrium free energy 
curves. Remarkable agreement is obtained between the 
two approaches. Also, far from negligible, _B C orr results 
in a B about 40% smaller than J3 v j r , as can also be gath- 
ered from TablelH The leading term to the correction, the 
change in /i c with density, suggests that deletion of lat- 
tice sites weakens the system's response to compression. 
The changes in lattice site occupancy permitted by parti- 
cle overlap thus increase the crystal compressibility com- 
pared to simple affine transformations. Generally, this 
still translates into an increase of B with density. Note 
also that the temperature dependence is rather weak. For 
T = 0.5 — 1.1, the curves in Fig.[5]appear to collapse onto 
a master function, which suggests that entropic effects 





FIG. 3: (Color online) Density fluctuations in multiple oc- 
cupancy crystals either stem from (left) fluctuations in the 
unit cell volume or from (right) rearrangements in particle 
distributions between lattice sites. 



have little impact in this regime. 

If sound waves have a period shorter than the time it 
takes for particles to redistribute between unit cells, then 
we can distinguish between two different sound modes 
in cluster solids. Density fluctuations stem either from 
changes to the unit cell volume v c at fixed cluster oc- 
cupancy n c or from fluctuations in n c at fixed v c , as 
schematized in Fig. [31 At constant T and V, the free 
energy density fluctuations Af = AF/V are then 

A/ = cn{p c An c ) 2 + c 2 2{n c Ap c ) 2 + ci2p c n c An c Ap c , (5) 



where the elastic constants 
1 f dp 



en 
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P\dp c J „ 



can be obtained by numerically differentiating the simu- 
lation results. A change of variables to sound-mode space 
x± = y 1 c\\p c An c ± y/C22n c Ap c , diagonalizes the expres- 
sion 



A/=(l+/?)4 + (l-/3)a2 



(6) 



with coupling constant j3 — Cyif ^c\\C22- The elastic 
coefficients and the coupling constant are presented for 
three different state points in Table [H Though these re- 
sults are insufficient to paint the full physical picture, a 
couple of comments are in order. First, the en term cor- 
responds to the "permeation" of particles and is expected 
to be heavily damped, while the C22 term is the equiva- 
lent of a longitudinal sound wave, which will propagate 
for long wavelengths. Second, for the temperature and 
density range under study the first constant and the cross 
term increase only very little with T and p, while the C22 
more than doubles. The increased density from one state 
point to the next is most certainly responsible for that, 
since higher temperatures would tend instead to facili- 
tate lattice spacing fluctuations for a constant repulsive 



4 



T P 




Cll C12 C22 P 


0.5 4.3 
0.8 6.2 
1.1 8.2 


48.2 89.8 41.6 

100.5 177.1 76.6 

176.6 308 131.0 


1.335 2.95 10.03 0.805 
1.346 3.05 15.35 0.670 
1.350 3.10 21.0 0.582 



TABLE I: Bulk modulus decomposition and the sound mode 
elastic and coupling constants for three different multiply- 
occupied crystal state points. 
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energy barrier. To the best of our knowledge, no theo- 
retical predictions exist with which to further compare 
these results. 

Starting from the formalism developed by Swope and 
Andersen [lj, we have presented how simulations and 
experiments of multiple-occupancy crystals critically de- 
pend on the chemical potential associated with the in- 
sertion of a lattice site. Taking this into account within 
a simulation allows for the precise determination of the 
equilibrium phase diagram of cluster crystals, which 
is much more subtle than for the traditional, single- 
occupancy sort. Also, even though the chemical poten- 
tial associated with lattice site insertion is strictly zero in 
equilibrium, its constrained derivatives are not. This has 
considerable impact on the calculation the bulk modulus 
and the two sound modes' elastic constants, for exam- 
ple. Departing more drastically from thermodynamics, 
long-lived non-equilibrium structures of cluster-crystal 
forming dendrimers might even be observable in rapidly 
quenched experimental systems, if the resulting ordered 
solids happen to end up in states with fi c ^ 0. These 
metastable crystals would then undergo phase transitions 
at different state points than those predicted by equilib- 
rium thermodynamics. Finally, the generalization of the 
free energy calculation methodology presented here has a 
broader applicability than for multiply-occupied crystals. 
It would also be the appropriate way to simulate systems 
with variable occupancy of lattice sites, such a micellar 
crystals or microphase separated colloids. 
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